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Abstract. We present a dynamical analysis of 673 galactic 
Planetary Nebulae, using a two-integral axisymmetric model 
with a Kuzmin-Kutuzov Stackel potential. The method fits the 
kinematics to the projected moments of a distribution func- 
tion, by means of Quadratic Programming. The 2.2 jjia COBE 
brightness map has been used after correction for the interstel- 
l/^ lar extinction as a projected star counts map in the modeling, 
0^ because it constitutes a galactic distribution view of evolved 
red populations which are considered to be the progenitors of 
^ PNe. 

^ The model we have obtained provides a 2-integral distribution 
0^ function for the COBE 2.2 jjia map, and thus a fortiori a de- 
projection of it, which allows moreover the identification of all 
the major Galactic components. We derive the density laws 

00 for them. The projected velocity dispersions are not well fit- 
ted though, especially in the disk, which points at the likely 
presence of a third integral. If this result can be confirmed by 
additional data, this would mean that for the first time the 

C ^ ^ presence and importance of a third integral on a global scale 

04 is demonstrated. 
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Q 1. Introduction 

H 

The dynamical modeling of various stellar populations in our 
Galaxy is an important tool in the study of its formation and 
evolution. Any such sample of course must consist of selected 
candidates, which often restrict our attention to luminous stars 
such as K- and M- Giants (Terndrup 1993), variable stars such 
as Mira's (Whitelock, 1993), post-AGB stars such as OH/IR 
stars (Sevenster, Dejonghe & Habing 1995) or Planetary Neb- 
ulae, hereafter PNe (AUer, 1993). None of these tracer objects 
necessarily belong to the main-stream pack, which complicates 
the interpretation of the data. 

Dynamical studies of the Galaxy as a whole need reasonably 
homogeneous samples over large volumes. At present, only the 
OH/IR stars (te Lintel Hekkert, 1990) and the PNe come into 
consideration for this. The OH/IR stars, which are strong in- 
frared emittors, oxygen-rich post-AGB stars in the final state 



Send offprint requests to: S. Durand 



of the red giant phase are especially suited, since the avail- 
able samples are products of blind surveys or result from well- 
defined selection procedures on otherwise blind surveys (no- 
tably IRAS), as demonstrated by Pottasch et al. (1988) and 
Ratag et al. (1990). As to the PNe, they are transient objects 
of considerable astrophysical interest. Their existence is based 
upon the synchronous evolution of their 2 components : a hot 
central star surrounded by a circumstellar nebula in expansion. 
PNe are thought to descend from low- and intermediate-mass 
stars (Mj < 8Mq) and probably cover a large range of ages. 
They represent a very short and very late phase in stellar evo- 
lution, located beyond the AGB stage. In their death agony 
they probably end up as white dwarfs. PNe can be tracked 
down relatively easily throughout the entire galaxy thanks to 
their rich emission line spectra that make them suitable as test 
particles for dynamical studies. 

2. The database 

Up to now, 1143 objects are classified as "true or probable" 
PNe in the Strasbourg- ESO Catalogue of Galactic PNe (Acker 
et al., 1992). Most of them seem to belong to a thick disk 
population, while, in accordance to the bulge selection criteria 
of the Catalogue {\l\ and |6| < 10°, radio fluxes F{6cm) < 
100 mJy and angular diameters < 20 arcsec), about 300 are 
probably bulge objects. 

We will use only the 673 PNe for which radial velocities are 
available. The database used in our modeling comprises two 
distinct data sets : 

2.1. Galactic positions and radial velocities of 673 PNe 

The radial velocities stem from miscellaneous spectroscopic 
surveys. Most are compiled originally in Schneider & Terzian 
(1983) and in the Strasbourg-ESO Catalogue (Acker et al., 
1992). Some others, about one hundred still non-published data 
with a precision of about 50 km/s, come from the bulge sur- 
vey of low spectral resolution performed by A. Acker and B. 
Stenholm (Acker et al., 1991). Moreover, about 60 radial ve- 
locities with a precision of nearly 1 km/s and originated from 
high resolution observations performed by A. Acker and A. Zi- 
jlstra were added to our sample. These last radial velocities 
constituted the 'reference set' for the calibration of the low- 
and medium-resolution ones (Durand et al., 1995). 
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Fig. 2. A comparison of PNe (dotted curves) witli tlie extinction corrected intensities provided by COBE at 2.2 /im (solid curves), the 
latter scaled so as to produce roughly the same numbers as the other two. Note the effects of the (optical) extinction in the PNe sample at 
small latitude, and the positive detection bias for Bulge PNe. 



About 300 of these 673 PNe are thought to belong to the galac- 
tic bulge. 

Since the Strasbourg-ESO catalog is to a large extent a com- 
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Fig. 1 . The galactic distribution (up) and the radial velocities versus 
galactic longitude (down) of our 673 PNe sample. 

pilation of all known or conjectured PNe in the literature, its 
selection biases are rather poorly understood. Moreover, the 
catalog draws mostly from optical work, and is thus severely 



incomplete in the plane. Therefore, it is not an easy-to-use 
instrument for producing star counts, which are affected by 
selection biases to zeroth-order and extinction effects. For ex- 
ample, we can easily identify a selection effect due to the special 
attention which the bulge attracts, causing a disproportional 
detection efficiency of PNe in that region (see also discussion of 
figure 2). We will thus use the catalog only to obtain the kine- 
matics of the PNe, i.e. projected mean streaming and projected 
velocity dispersion. 

2.2. COBE intensity map of the full sky at 2.2 jjm 

At 2.2 jUm, emissions of late-K and M giants dominate the spec- 
trum (Arendt et al., 1994). Our galaxy is furthermore relatively 
transparent to (cold) interstellar material in that K band, and 
thermal emission from dust is still negligible. Hence, the 2.2 
jUm COBE picture is a rather reliable map of the projected 
spatial distribution of galactic AGB stars. Moreover, there is 
at the present day a general "consensus that AGB stars are im- 
mediate progenitors of PNe" (Habing et al., 1993): red giants 
and PNe have evolutionary links through the OH/IR phase, 
although the precise nature of all the transition objects is not 
yet really known (Habing, 1990). So it is not unreasonable at 
first sight to use the 2.2 jjia COBE map as a substitute for 
the spatial distributions of galactic AGB and post- AGB popu- 
lations. We will now show that corroborating evidence for the 
validity of this procedure can be found in the comparison of 
the data sets. 

Clearly the raw COBE map must be corrected for interstellar 
extinction. Following Arendt et al. (1994), we used the relations 
between near IR colours (caused by extinction), as determined 
by Rieke & Lebofsky (1985). We first constructed a colour ex- 
cess map E{ J — K) by the combination of the two first near 
IR brightness COBE maps, and by choosing a unique intrin- 
sic colour ( J — K) for the Galaxy, thus neglecting the possible 
presence of a variation in the intrinsic colour between the in- 
ner and the outer parts of the Galaxy. We have then built an 
optical depth map of the galaxy at 2.2 jjia by associating the 



reddening law. Finally we transformed the 2.2 jjia brightness 
COBE map into a counts density map by a simple scaling, so 
as to reproduce roughly the same numbers as the stars sam- 
ple. Figure 2 compares the PNe counts and the scaled 2.2 jjia 
COBF map. The first row of plots clearly shows the effect of 
the optical extinction, which affects the PNe counts. At larger 
latitudes however, the curves are very similar. At the present 
day it's hardly possible to correct for the extinction in these 
stars counts. The second row reveals that, even in the pres- 
ence of extinction at 6 = 2°, the slopes remain very similar for 
\l\ > 20°. On the other hand, our PNe counts appear effectively 
overestimated in the bulge, and that can be well understood 
since we have added to the original PNe catalog nearly a hun- 
dred of new bulge measurements. 

From this we see that it is certainly reasonable to use the 
COBF map as a substitute for the PNe counts, in order to 
overcome the different biases in the PNe catalogs. This can be 
asserted, not only for general reasons of astrophysical nature, 
but also on the basis of the comparison of both data sets, there 
where it makes sense to compare them. 

We have restricted our study to the galactic region 2 < |6| < 8 
for the following reasons: 

- |6| < 2 is likely to be a very complicated region where opti- 
cal depths reach their highest values, in so far as the extinc- 
tion corrected COBF map cannot reasonably represent the true 
galactic K emission. 

- PNe data at |6| > 8 are unfortunately not numerous enough 
to allow a correct estimation of the different moments of the 
velocity distribution. 

Moreover, the contribution of the zodiacal light is minimal in 
the K band as evidenced in the composite 1.25, 2.2 and 3.5 
jUm COBF images (Hauser 1993), so we have neglected it as 
is common practice in near-infrared studies of the Galactic 
emissivity distribution (see also Kent et al., 1991). In our case 
in particular, this approximation is not critical, since the effect 
reaches levels that are below the accuracy we can reasonably 
obtain for the 0-th moment of the distribution function. 



3. The modeling 

3.1. The framework 

PNe are highly evolved objects supposed to represent a dynam- 
icaly relaxed population. Therefore, they are in a sufficiently 
steady state to legitimate the application of an equilibrium 
model. One of the first choices to be made is the modeling ge- 
ometry: we chose here the axisymmetric approximation. This 
could appear at first sight as rather outdated, because of the 
mounting evidence that the Galaxy is a complex triaxial sys- 
tem (e.g. Spergel D.N., 1992). However, the triaxiality in PNe 
kinematics is not obvious, if not absent, so it seems legitimate 
to us (by Occam's razor) to firstly exploit the axisymmetric 
geometry. 

Our basic aim is the construction of the distribution function 
of the PNe population, which is simply the probability density 
of the sample in the 6 dimensional phase space (x, v). It will be 
later of great help to appreciate the different orbital structures 
present in the PNe population (section 4). 

The potential we use in this work is a Kuzmin- Kutuzov Stackel 
potential with a halo-disk structure (Batsleer & Dejonghe 
1994). Its general expression in spheroidal coordinates can be 
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d and h standing for 'disk' and 'halo' 

A, spheroidal coordinates 

G the constant of gravity 

M the total mass of the galaxy 

K = Md/M 

q a parameter defined by g = Ad — 



i^d 



iyh>0. 



Let's recall briefiy that spheroidal coordinates (A, i^) are de- 
fined as the roots for r of the equation : 
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where (r, (j), z) are the cylindrical coordinates and a, j negative 
constants. The axis ratios - of each of the two coordinate sur- 

c 

faces (with a = —\fa and c = — t/v) and the K parameter are 
selected in order to reproduce the shape of the galactic rota- 
tion curve with its characteristic fiat extent. For this, we use 
for our mass model: 
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The selection of the parameters is based upon the requirements 
that the rotation curve be fiat, that the radial and thickness 
scalelengths of the disk are close to the commonly accepted 
values, and that the surface and spatial densities of the disk in 
the solar neighbourhood is reproduced (for more details about 
the selection of the parameters, see Batsleer & Dejonghe, 1994). 
Note that the shape of the halo is almost spherical and that the 
adopted K parameter implies a disk mass of 2.94 X 10^" Mq. 
The (stationary) distribution function is usually described via 
Jeans' Theorem as a function of integrals of motion which 
are, in the present case, the orbital binding energy E, the z- 
component of the angular momentum Lz, and a third integral 
which is simple and analytic in the Stackel case. Since informa- 
tion about vertical motions is still practically non-existent, the 
explicit presence of a third integral in the distribution function 
will not directly reveal itself in the kinematics. Therefore, only 
the failure of a complete dynamical 2-integral model to fit the 
data may point to the presence of the third integral. The ma- 
jor approximation made with the utilization of a two-integral 
model is that radial and vertical dispersions are supposed to 
be equal; this can be easily checked just by looking at the ex- 
pressions, in cylindrical coordinates, of the two first integrals 
of motion E = !/'(r, z) — \{vl- + f e + *'z) Lz = r vg. The 
contributions of the Vr and Vz components are completely iden- 
tical contrary to the vg component which appears in an inde- 
pendent way in Lz- That 2-integral approximation should be 
rather good in the thin disk where dispersions are thought to 
be lower than anywhere else since orbits there are nearly circu- 
lar, but problems should arise elsewhere, and more particularly 
in the thick disk. 




Fig. 3. A comparison of tlie data (solid curves) and of our best 2 integral model (dashed curves) at different galactic positions. 

From left to right : the projected counts per square degree, the projected mean velocities (in km/s) and the projected velocity dispersions 

(in km/s). The error bars on the PNe kinematics have been obtained by varying the sampling strategy. 



3.2. The QP method 

In this method, the distribution function is written as a sum 
of elementary distribution functions, 

F{E,L,) = '^c,F,{E,L,) with I = 1(0,13,-1), (1) 

which are of 2 kinds (families) that we call ^building blocks^: 

1. F^AE,L) = E^Ei^Y, 

which are essentially Fricke components associated with fairly 
slowly rotating models with high velocity dispersion (halo- 
type). 

2. Ec,i3,-y{E, L) = exp{-^) S"^ {2SL^ f {E - SJ' if S < E 
E„},l,(E,L)=0 ifS>E 

with S the minimum orbital energy of an orbit which can bring 
a star with an associated angular momentum Lz at a given 
height zo. 

The second building block produces models with much smaller 
scaleheights than the first one, and with nearly circular orbits 
(disk-type). Note that components with /3 > have a zero 
central density, and that the presence of the exponential factor 
is useful to control the radial behaviour of the mass distribution 
in the disk (Batsleer & Dejonghe, 1995). 

In order to determine the coefficients c, of the distribution 
function (1), we consider a grid in galactic coordinates on which 
data points are placed: each of the grid points is associated 
with the O"*, 1*' and 2"'' moments of the projected velocity 
distribution, calculated by taking into account the 15 stars 
closest to the grid point (see figure 4 and table 1). 
This yields for each point on the sky 3 observables iJobs{x,v) 
that can be estimated by means of a function ij{x,v) : 



IJ.mod{xi,vi) = J ij{x,v) E{I)dx dv 

= ^^c, / Hmod(x, v) Ei(I) dx dv = CtlJt{xi, Vl). 
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Fig. 4. The grid used for the evaluation of the different moments 
of the projected velocity distribution, superposed on nearly all our 
data points, with positions that were flipped into the first galactic 
quadrant. 



In our case r; is the product of the full velocity space available 
at a spatial point and the line of sight. Then we can construct 



s X "type function 



Wl[Hobs(xi, Vl) - ^C,jU,(2;i, vi)f 



with wi weights that can be freely chosen. This function is min- 
imized during the calculations with additionnal constraints to 
force a positive distribution function in all phase space (De- 
jonghe, 1989). This is a Quadratic Programming problem. This 
method has the advantage to accomodate inhomogeneous data 
sets (photometry, counts, kinematics) and additionnaly im- 
posed constraints, as long as these can be put in a form which 
is linear in the distribution function. Further details of the 
method as applied on similar problems can be found in Bat- 



Table 1. counts (in 10 ^ counts/sq degree) providing by COBE, projected mean velocities (in km/s) and pro- 
jected mean dispersions (in km/s) coming from PNe data, all given for each grid points (l,b). 
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sleer & Dejonghe (1995) and Sevenster, Dejonghe & Habing 
(1995). 

4. Results on PNe dynamics 
4.1. The validity of the modeling 

The quality of the fit can be checked directly by comparing the 
data and the model, as shown in figure 3 for different galactic 
regions. We see that the model fits perfectly the COBE density 
and the projected mean velocities. On the other hand, there is 
an obvious discrepancy in the fit of the projected dispersions, 
particularly in the disk where the errors are nearly of the same 
scale as the data : this is probably the signature of an explicit 
third integral dependence in the distribution function (see 3.1). 

4-2. The orbital structures 

4.2.1. The distribution function in integral space 

A schematic representation is shown in figure 5 and the PNe 
application in figure 8 (top right). Both figures are shown in 
the coordinate space (Lz\/2E, E), and not (Lz, E) in order to 
keep integral space finite. For a given point {r,z), the region 
of the allowed orbits is bounded by an ellipse which is tangen- 
tial in two points at the curved outer boundaries which are 
the loci of the circular orbits (one point for prograde rotating 
populations, at Lz > 0, and another for retrograde rotating 
populations, at Lz < 0). We recall that the positive Lz axis 
points towards the galactic north pole by convention, so our 
Galaxy is rotating retrograde. Each circular orbit curve (which 
occurs by definition at z=0) is shown as far as r = 15 kpc. 

4.2.2. The distribution function in turning point space 

This representation has the advantage that it is more intuitive 
than the previous one. A schematic view is shown in figure 
6 and the PNe application in figure 8 (top left). The coordi- 
nate system is (rp,ro), being respectively the pericenter and 
the apocenter of the orbit. The loci of the circular orbits are 
defined by the hues rp = -\-r a (prograde orbits) and fp — — fa 
(retrograde orbits), and the allowed orbits are thereby placed 
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Fig. 5. A schematic representation of the distribution function in 
integral space (L^V 2-£/, E) displayed up to a distance r of 15 kpc 
from the galactic centre. The two curves E = Sq[Lz) correspond to 
the circular orbits (occuring by definition at ^ = 0). The two curves 
E = Sz^{Lz) correspond to the minimum binding energy that can 
bring a star at most as high as zq . 



in the regions where > |rp| ; by convention rp has the sign 
of Lz. 

The 2 plates at the top of figure 8, presented in logarithmic 
density scale, show two very beautiful distribution functions 
which are quasi dominated by the second building-block family 
(disk-type). We note that most of the orbits are nearly circu- 
lar in the disk, that the bulge is not excessively populated (the 
orbits are very dispersed), and that the phase space density of 
the old disk increases from r = 5 kpc to about r = 10 kpc, 
whereafter it decreases smoothly. 
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Fig. 7. Cuts of the deprojected COBE map in some galactic regions. Tlie derived scale parameters are indicated. 
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Fig. 6. The distribution function in 'turning point' space 



4.3. The COBE map deprojection 

Our dynamical modeling also provides a method for depro- 
jecting an observed light distribution. It suffices to weigh only 
lightly the kinematic data in the x variable. The main ad- 
vantage of such a procedure is that we then know that the 



deprojected light distribution can be produced by orbits, i.e., 
stems from a positive distribution function. The deprojection, 
shown in figure 8 (bottom), assumes a galactocentric distance 
of 8 kpc. The bulge and old disk shapes are very clearly present. 

First, let's consider the bulge, the density of which is mainly 
concentrated inside 0.8 kpc. The radial profile is exponential 
(except for very small radii, of course), with a scale length Hr = 
0.4 kpc. The determination of the vertical scale height is more 
problematic because of the uncertainty in the correction of the 
interstellar extinction : it varies smoothly from = 180 pc at 
Z = kpc to = 270 pc at Z = 0.8 kpc (see figure 7). If we 
take 235 pc as the average value of H^, then the corresponding 
bulge scale height ratio is 0.56, in very good agreement with 
other axisymmetric bulge axis ratios values (0.6 in Weiland et 
al. (1994), who use the ratio of the average projected 2.2 jum 
COBE intensities, or 0.53 in Dwek et al. (1995), who use a 
Gaussian-type fit on the observed 2.2 jum COBE map, or again 
0.61 in Kent et al. (1991), using their oblate spheroid model), 
but lower in general than the ratios in triaxial bulge models. 

The bulge outer region at large ^-height, which probably forms 
a smooth transition towards the halo, clearly follows a power 
law of index 3.2 (figure 7). 

As to the disk, the resulting exponential scale length is 2.6 kpc 
(at the bottom of figure 8). Here there seems to be an obvious 
chromatic dependence: 2-3 kpc in near IR, 3.5 - 5.5 kpc in 
optical and 4.5 - 6 kpc for IRAS OH/IR stars (Kent et al., 
1991). Kinematical determinations show also a large spread 
but generaly favour a larger scale length (4.4 kpc for Lewis 
& Ereeman, 1989 ; 4.3 - 5.8 kpc for van der Kruit, 1989). On 
the other hand, our small value is in agreement with recent 
determinations suggesting that the disk scale length is indeed 



Fig. 8. At the top: logarithmic-scaled images of the distribution function in integral space (right) and in turning point space (left). Black 
indicates lower intensity and white indicates higher intensity. 

At the bottom: The 2-integral deprojection of COBE map, presented in logarithmic density scale. 



rather short: Robin et al., (1992), have found Hf = 2.5 kpc 
using UBV photometry toward the anticentre down to m = 25 
; Fux & Martinet, (1994), have resolved the asymmetric drift 
equation by substituting the velocity dispersion gradients with 
analytical approximations and an application to old disk kine- 
matical data in the solar neighbourhood (with a constant Hf ) 
has lead them to a value of 2.5 kpc. 

The exponential scale height Hf given by the model is constant 
and equal to 150 pc up to r = 5.5 kpc, and increases monoton- 
ically beyond that radius, with a gradient of about 20 pc/kpc, 
to attain for example 210 pc in the solar neighbourhood. This 
gradient is quite similar to the value found originally by Kent et 



al. (1991), and its presence has been recently suggested again 
(Fux & Martinet, 1994 ; Ibata & Gilmore, 1995). As to the 
small value of the scale height, it may point to the fact that 
the 2.2 jUm emission detected by COBE does not solely origi- 
nate from an old population, but that a considerable fraction 
of it actually stems from young disk populations. 
There is another (non negligible) disk-like component in our 
configuration space density map at ^ > 0.3 kpc. It presents 
all the characteristics of a thick disk, which was originally pro- 
posed by Gilmore et al., (1983), and is thought to be populated 
by intermediate populations. Besides, we have attempted to 
make a global estimation of the averaged ages of the disks 



vided by Widen (1977). The corresponding 2-integral depro- 
jection of the dispersions is shown in figure 9, and show up the 
possibiUty that PNe in the thick disk may be estimated to be 
around 7-9 Gyrs old, contrary to the young disk ones where 
averaged ages may turn around 2-5 Gyrs. 

The radial scale length of the thick disk here appears to be 
consistent with that of the old disk, namely 2 kpc. Vertical 
scaleheights are many, with typical values varying between 0.5 
and 1 kpc (figure 7). 

Finally, at about z = A kpc a transition between the thick disk 
and halo occurs, which is probably not well determined by our 
2-integral model. 
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Fig. 9. Radial (up figure) and tangential (down figure) dispersions 
of PNe at Z = 0.2, 1, 2 and 3 kpc, as provided by the 2-integral 
deprojection. 



distribution function is made of a PNe sample. 

An important conclusion of our paper remains that a 2-integral 
model seems not really adequate to characterize the dynamical 
state of PNe, thereby demonstrating the possible presence and 
importance of the third integral on a global scale. Of course, 
the dependence of the distribution function on the third in- 
tegral in the solar neigbourhood has been established a long 
time ago, and it is widely believed on firm theoretical grounds 
that this is also the case on a global scale. The observational 
demonstration is not trivial however, since the effect of the 
third integral is mostly felt in the vertical velocity, which con- 
tributes only weakly to the radial velocities of objects in the 
disk due the position of the sun. Global models of the kind 
presented in this paper are able to exclude 2-integral models, 
because for such models Ur = <Tz, fr is well observable in the 
radial velocities and both dispersions are uniquely determined 
by the morphology of the stellar component (i.e. the COBE 
map). 

The failure of our 2-integral models is also confirmed by the 
results of Sevenster et al. (1995), where the same modeling 
scheme is applied on a sample of 700 galactic OH/IR stars (the 
sample 'tLH'). If we compare the lower right panel of their fig- 
ure 4 with the lower right panel of our figure 3, we see a very 
similar behaviour, though probably a bit less pronounced for 
the OH/lRs: the velocity dispersion of the models are too low 
at large longitudes. Sevenster et al. also note in passing that 
their tLH sample may need an additional integral. This is cer- 
tainly to be expected, since the kinematics of the OH/lRs and 
the PNe are very similar, except that the dispersions in the 
OH/IR population are a bit lower and that the modeling of 
the OH/lRs used star counts instead of the COBE map be- 
cause the OH/IR sample is much more homogeneous. 
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5. Concluding remarks 

The extinction-corrected 2.2 jum COBE map proves to be a 
very useful substitute for AGB and post-AGB star counts in 
the presence of non-kinematic detection biases. It can also be 
used to isolate in the models the effects of kinematics alone, 
and one could thus compare the dynamics of different samples 
if one uses systematically the same projected density map for 
all of them. 

The bulge-disk decomposition is very obvious in the density 
map in (r, z) obtained via the 2-integral deprojection of the 
COBE data, which allows us to estimate the scale parame- 
ters of all the principal galactic components: for example, our 
model conforts some recent ideas about the disk as a scale- 
length rather short (2.6 kpc) and a vertical profile different 
from 'pure-disk' galaxies ones where scale heights are indepen- 
dent of galactocentric distance. There is also evidence that the 
2.2 jUm disk emission has a considerable component that orig- 
inates from young populations. 
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